# MC simulation ambiguity model

rm(list = ls())
library(foreign)

# utility functions and voting probabilities
uij <- function(xi,zj,sigmaj) {
  if (length(xi) > 1 | length(zj) > 1) {
    n.x <- length(xi)
    n.z <- length(zj)
    x <- matrix(rep(xi,n.z),n.x,n.z,byrow=FALSE) 
    z <- matrix(rep(zj,n.x),n.x,n.z,byrow=TRUE)
    sigma <- matrix(rep(sigmaj,n.x),n.x,n.z,byrow=TRUE)
    t <- -sigma^2 - (x-z)^2
  } 
  else {
    t <- -sigmaj^2 - (xi-zj)^2
  }
  t
} 
rohij <- function(xi,zj,sigmaj,z,sigma) {
  exp(uij(xi,zj,sigmaj)) /  rowSums(exp(uij(xi,z,sigma)))
}
roh <- function(x,z,sigma) {
  exp(uij(x,z,sigma)) /  rowSums(exp(uij(x,z,sigma)))
}
u2j <- function(xmj,hetj,zj,sigmaj) {    
  if (zj <= xmj) {
    lower <- xmj-sqrt(3)*hetj
    upper <- zj + sqrt(3)*sigmaj
    mu.post <- (upper+lower)/2
    sigma.post <- (upper-lower) / (2 * sqrt(3))
  }
  if (xmj < zj) {
    lower <- zj-sqrt(3)*sigmaj
    upper <- xmj + sqrt(3)*hetj
    mu.post <- (upper+lower)/2
    sigma.post <- (upper-lower) / (2 * sqrt(3))
  }
  if (lower > upper) {
    u <- NA
  }
  else {
    u <- -sigma.post^2 - (xmj-mu.post)^2   
  }
  u
}
roh2j <- function(xmj,hetj,zj,sigmaj) {
  exp(u2j(xmj,hetj,zj,sigmaj)) / ( exp(u2j(xmj,hetj,zj,sigmaj)) + exp(uc) )  
}
p2 <- function(z,sigma) {
  t <- rep(0,length(z))
  for (j in 1:length(z)) t[j] <- roh2j(xm[j],hetj[j],z[j],sigma[j])  
  t
}
# contraction mapping T
g <- function(x,z,sigma,EM) {
  t <- (EM <= z)*2-1
  zstar <- (  6 * colSums( (roh(x,z,sigma)-roh(x,z,sigma)^2) * x) + 
       alpha * length(x) * (p2(z,sigma)-p2(z,sigma)^2) *
       ( (-t) * 2* sqrt(3) * sigma + t * sqrt(3) * hetj + 2* xm   )    ) /
    (   2 * alpha * length(x) * (p2(z,sigma)-p2(z,sigma)^2) + 6* colSums(roh(x,z,sigma)-roh(x,z,sigma)^2)  )
list("z"=zstar)
}
h <- function(x,z,sigma,EM) {
  sstar <- (   alpha * length(x) * (p2(z,sigma)-p2(z,sigma)^2) * 
       (  3 * hetj + 2 * sqrt(3) * abs(xm -z))  ) /
    (   6 * (alpha * length(x) * (p2(z,sigma)-p2(z,sigma)^2) + colSums(roh(x,z,sigma)-roh(x,z,sigma)^2)  ) ) 
list("sigma"=sstar)
}
# Nash equilibrium
ne <- function(x,start.z,start.sigma,precision) {
  i <- 1
  z <- start.z
  sigma <- start.sigma
  d.z <- 1000
  d.sigma <- 1000
  h.z <- 0
  h.sigma <- 0
  while ((d.z > precision | d.sigma > precision) & i < 500) {
    i <- i+1
    z_1 <- z
    sigma_1 <- sigma
    okay <- FALSE
    ii <- 0
    while (okay == FALSE & ii < 100) {
      ii <- ii + 1
      EM <- colSums( (roh(x,z_1,sigma_1)-roh(x,z_1,sigma_1)^2) * x) / colSums((roh(x,z_1,sigma_1)-roh(x,z_1,sigma_1)^2))
      z <- g(x,z_1,sigma_1,EM)$z
      okay <- all(!is.na(z))
      if (okay==FALSE) {
        sigma_1[is.na(z)] <- sigma_1[is.na(z)]*1.1
      }
    }
    hh <- h(x,z_1,sigma_1,EM)
    sigma <- hh$sigma
    d.z <- sum(abs(z-z_1))
    d.sigma <- sum(abs(sigma-sigma_1))
    h.z <- c(h.z,d.z) 
    h.sigma <- c(h.sigma,d.sigma)
  }
  list("z"=z,"sigma"=sigma,"it"=i,"start.z"=start.z,"start.sigma"=start.sigma,"EM"=EM)
}

# simulation with 100 runs, 10 starting values
# number of simulations: sims
sims <- 100
n.start <- 10
# fixed parameters
uc <- -1
n <- 1000
alpha <- 1
prec.pow <- 5
prec <- 10^(-prec.pow)
# matrices to collect results
mc <- array(list(NULL), c(sims,10))
mc.xm <- mc.x <- array(list(NULL), c(sims,1))
# uni.z.ne <- NULL
# uni.sigma.ne <- NULL
# uni.count.ne <- NULL
# Monte Carlo simulation, sims runs
for (ii in 1:sims) {
  cat(".")
  x <- rnorm(n,0,1)
  p <- sum(rmultinom(1,size=1,prob=c(1,1,1)) * c(1,2,3))*2 + 1
  distinct.xm <- p-1
  while (distinct.xm < p) {
    xm <- round(rnorm(p,0,1),digits=3)
    distinct.xm <- length(table(xm))
  }
  xm <- sort(xm)
  hetj <- rep(0.2,p)
  z.ne <- NULL
  sigma.ne <- NULL
#   it.ne <- NULL
#   start.z.ne <- NULL
#   start.sigma.ne <- NULL
#   xm.ne <- NULL
#   xmean.ne <- NULL
  # run for n.start different starting values
  for (i in 1:n.start) {
    start.z <- runif(p,0,1)*xm
    start.sigma <- (1/sqrt(3)*abs(start.z-xm)+1/2*hetj)*runif(p,0,1) 
    t1 <- ne(x,start.z,start.sigma,prec)
    z.ne <- rbind(z.ne,round(t1$z, digits=prec.pow))
    sigma.ne <- rbind(sigma.ne,round(t1$sigma, digits=prec.pow))
#     it.ne <- rbind(it.ne,round(t1$it, digits=prec.pow))
#     start.z.ne <- rbind(start.z.ne,round(start.z, digits=prec.pow))
#     start.sigma.ne <- rbind(start.sigma.ne,round(start.sigma, digits=prec.pow))
#     xm.ne <- rbind(xm.ne,xm)
#     xmean.ne <- rbind(xmean.ne,mean(x))
    mc[[ii,i]] <- t1
  }
  mc.xm[[ii]] <- list("xm"=xm)
  mc.x[[ii]] <- list("x"=x)
  # plot simulation for different starting values  
  par(mfrow=c(2,1)) 
  par(mar=c(2,1,2,1))
  par(oma=c(1,1,1,1))
  dotchart(z.ne, labels="",main="z",xlim=c(-3,3))
  dotchart(sigma.ne,labels="",ylab="sigma",main="sigma",xlim=c(0,3))
  # collect
#   uz <- unique(round(z.ne,prec.pow-2)) 
#   us <- unique(round(sigma.ne,prec.pow-2))
#   uz <- cbind(uz,matrix(0,nrow(uz),11-p))
#   us <- cbind(us,matrix(0,nrow(us),11-p))
#   uzz <-cbind(seq(1,nrow(uz)),uz)
#   uni.z.ne <- rbind(uni.z.ne,cbind(rep(ii,nrow(uzz)),uzz)) 
#   uss <-  cbind(seq(1,nrow(us)),us)
#   uni.sigma.ne <- rbind(uni.sigma.ne,cbind(rep(ii,nrow(uss)),uss))
#   uni.count.ne <- rbind(uni.count.ne,nrow(uz))
}

# save.image(file = "simdata.RData")
# load(file = "simdata.RData")
# 
# # figure for paper
# cases <- sample(c(1:sims),4)
# par(mfrow=c(2,2)) 
# for (t in cases) {
#   q <- mc[[t,1]]
#   qq <- mc.xm[[t]]  
#   qqq <- mc.x[[t]]
#   PP <- length(q$z)
#   p <- seq(1,PP)
#   rgx <- c(-2,2)
#   rgy <- range(c(0,PP))
#   par(mar=c(2,4,1,1))
#   plot(rgx,rgy, type="n", axes=F, xlab="", ylab="", cex.lab=1,font.lab=1 )
#   par(new=TRUE)
#   plot(qq$xm,p, xlim=rgx, ylim=rgy, pch=20, xlab="",ylab="",yaxt="n")
#   for (ii in 1:PP) {
#     lines(cbind(qq$xm[ii]-sqrt(3)*hetj[1],qq$xm[ii]+sqrt(3)*hetj[1]),c(ii,ii)) 
#   }
#   par(new=TRUE)
#   plot(q$z,p-.15, xlim=rgx, ylim=rgy, pch=8, xlab="",ylab="",axes=FALSE)
#   for (ii in 1:PP) {
#     lines(cbind(q$z[ii]-sqrt(3)*q$sigma[ii],q$z[ii]+sqrt(3)*q$sigma[ii]),c(ii-.15,ii-.15)) 
#   }
#   par(new=TRUE)
#   plot(q$EM,p, xlim=rgx, ylim=rgy, pch=1, xlab="",ylab="",axes=FALSE)
#   par(new=TRUE)
#   plot(density(qqq$x),axes=FALSE,main="",ylab="",xlim=rgx)
#   IV <- abs(q$EM-q$z)*sqrt(3)
#   reg <- round(lm(q$sigma ~ IV)$coef[2],2)
# }


